Leaflet Homework Key

library(dplyr)
library(sf)
library(tidyverse)
library(leaflet)
library(sf)
library(leaflet.providers)
library(RColorBrewer)
library(tigris)

Dot Density Plot

# education_url <- "https://github.com/dyerlab/ENVS-Lectures/raw/master/data/edu_attainment-shp.zip"
# 
# file <- tempfile(pattern = "edu", tmpdir = ".", fileext = ".zip")
# file
# 
# download.file(education_url, destfile = file)
# unzip(file, exdir = ".")
# 
# 
# points <- st_read("edu_attainment-shp")

setwd("~/Spring2021/GISinR/GISinR")

The working directory was changed to C:/Users/Charis/Documents/Spring2021/GISinR/GISinR inside a notebook chunk. The working directory will be reset when the chunk is finished running. Use the knitr root.dir option in the setup chunk to change the working directory for notebook chunks.

points <- st_read("edu_attainment-shp")

Reading layer `edu_attainment' from data source `C:\Users\Charis\Documents\Spring2021\GISinR\GISinR\edu_attainment-shp' using driver `ESRI Shapefile'
Simple feature collection with 15745 features and 1 field
geometry type:  POINT
dimension:      XY
bbox:           xmin: -77.59891 ymin: 37.44779 xmax: -77.38568 ymax: 37.60228
geographic CRS: WGS 84

Convert education levels to factor - this allows us to use colorFactor() to create our color palette

points$education <- factor(points$education, levels = c("No HS diploma", "HS, no Bachelors", "Bachelors" , "Post-Bachelors"  )) 

Create color palette

factpal <- colorFactor(brewer.pal(n = 4, name = "RdBu"), points$education) #used to color categorical data

Create dot density plot

leaflet(points) %>% 
  addProviderTiles(providers$CartoDB.Positron)%>% #add provider tiles
  addCircles(stroke = FALSE,  fillOpacity = 1, #add points
              color = ~factpal(education)) %>% #color based on educational attainment
  addLegend(pal = factpal, values = ~education, title = "Education Attainment")%>% #add legend
  addScaleBar(position = "bottomright") #add scale bar

Choropleth plot

load in the census tract data

tracts <- tracts("VA", "Richmond city") %>%
  st_transform(4326)

Using FIPS code '51' for state 'VA'
Using FIPS code '760' for 'Richmond city'

Data Prep

edu_tracts <- tracts %>%
  st_buffer(dist = 0) %>% # buffer by 0 to avoid error message
  st_intersection(points) # intersect the point data with the tract polygons

edu_tracts %>%
  group_by(NAME, education) %>% # group by tract name and education level
  tally() -> tally # tally() counts the number of rows for the groups, so output is one row for every attainment level for each tract (so if a tract has all 4 attainment levels, there are 4 rows of data for that tract)

tally <- as.data.frame(tally) # converted to data.frame to be able to drop the geometry because the `pivot_wider()` function would not work the way I wanted it to when the geometry was present

tally <- tally %>%
  select(NAME, education, n)%>% # do not want the geometry 
  pivot_wider(names_from = education, values_from = n) # pivot the table so that instead of having multiple rows per tract, you have one row per tract with a column for each attainment level

tally[is.na(tally)] <- 0 # define NA values as zero

tally

Create attainment metric for choropleth

tally <- tally %>%
  mutate(edu = (((`No HS diploma` * 1) + (`HS, no Bachelors` * 2) + (`Bachelors` * 3) + (`Post-Bachelors` * 4))/
                  (`No HS diploma`  + `HS, no Bachelors` + `Bachelors` + `Post-Bachelors`)))

# I created a column, edu, and did a weighted average for each row. So 'No HS diploma' was weighted as 1, 'HS, no Bachelors' as 2, ect. 

# I used this new column as my average education attainment value to color my polygons

tally %>%
  select(NAME, edu) -> tally # no longer need the individual count data

Join the education data back with the tract data

tracts %>%
  left_join(tally, by = "NAME") -> tracts

Create the palette

pal <- colorNumeric(
  palette = "RdBu",
  domain = tracts$edu) # This time I want a continuos palatte, so I used `colorNumeric()`

Make the choropleth map

leaflet(tracts) %>%
  addProviderTiles(providers$CartoDB.Positron)%>% # add provider tiles
  addPolygons( fillOpacity = 1, # add tract polygons
               fillColor = ~pal(edu), smoothFactor = 0.2, # fill color based on the calculated 'edu' metric
               color = "grey", 
               opacity = 1,
               weight = 1) %>%
  addLegend(pal = pal, values = ~edu, title = "Education Attainment") # add legend

I did not like this legend, as it’s unclear what the numbers mean. With some googling I found a solution

labels <-c("No HS diploma (1)", "HS, no Bachelors (2)",
           "Bachelors (3)" , "Post-Bachelors (4)"  ) # define the labels for the legend

leaflet(tracts) %>%
  addProviderTiles(providers$CartoDB.Positron)%>% # add provider tiles
  addPolygons( fillOpacity = 1, # add tract polygons
               fillColor = ~pal(edu), smoothFactor = 0.2, # fill color based on the calculated 'edu' metric
               color = "grey", 
               opacity = 1,
               weight = 1) %>%
  addLegend(pal = pal, values = ~edu, 
            opacity = 1,
            title = "Education Attainment",
            labFormat = function(type, cuts, p) {  # Here's the trick
              paste0(labels)
            })

That looks better, I think

Combine the maps

Final map